%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Berechnung der optischen Blochgleichungen und Darstellung der Fluoreszenz,
%  das Ion ist mikrobewegt,
%  n = [r11 r22 r12 r21] ,
% --- |2> 
%  |
% --- |1> 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
I=1i/2;

G12=           1;															% decays P1/2 - S1/2
g12=         G12;															% decays of the coherences
o12=      .5*g12;															% Rabi frequencies
d12=           0;															% detunings
L12=g12/2+1i*d12;
mod=[1 3 10 30]*.1;

figure(1);					clf;
for kk=1:length(mod)
    m  =       mod(kk);														% modulation excursion
    f  =	        .5;														% modulation frequency
    D12=[-1:.01:1]*G12*10;

    A=[-G12		-I*o12	I*o12
       -1i*o12	-L12	   0
        1i*o12		0	-L12'  	];

    B=[0  		0		0
       0      1i*m		0
       0		0	  -1i*m		];

    b  =[0	-I*o12	I*o12		]';
    E  =1i*f*eye(3);

    nt=zeros(size(D12));
    for k=1:length(D12);
        d12=     D12(k);
        L12=g12/2+1i*d12;			A(2,2)=-L12;			A(3,3)=-L12;

        Sp=eye(3);					Sm=eye(3);
        for kk=15:-1:1																		%Generation of the continued fractions
            Sp=-(A+kk*E+B*Sp)\B;	Sm=-(A-kk*E+B*Sm)\B;
        end
        rst=-(A+B*(Sp+Sm))\b;
        nt(k)=rst(1);
    end
    hold on
    
    subplot(221),plot(D12/G12,real(nt));
    set(gca,'FontSize',16,'LineWidth',1);
    xlabel('\fontsize{16}\Delta_{12}')
    ylabel('\fontsize{16}{\it\rho}_{22}')
    % title(['\fontsize{16}\Omega/\Gamma=',num2str(o12/G12),'   Temperature='])
pause(.1)
end
hold off